Exploring tidal effects of coalescing binary neutron stars in numerical relativity 
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We study gravitational waves emitted in the late inspiral stage of binary neutron stars by analyzing 
the waveform obtained in numerical-relativity simulations. For deriving the physical gravitational 
waveforms from the numerical results, the resolution extrapolation plays an essential role for our 
simulations. The extrapolated gravitational-wave phases are compared with those calculated in 
the post-Newtonian (PN) and effective-one-body (BOB) formalisms including corrections of tidal 
effects. We show that the extrapolated gravitational-wave phases in numerical relativity agree well 
with those by the PN and EOB calculations for most of the inspiral stage except for a tidally- 
dominated, final inspiral stage, in which the PN and EOB results underestimate the tidal effects. 
Nevertheless, the accumulated phase difference between our extrapolated results and the results by 
the PN/EOB calculations is at most 1-3 radian in the last 15 cycles. 

PACS numbers: 04.25.D-, 04.30.Db 



I. INTRODUCTION 

The inspiral and merger of coalescing compact bina- 
ries are among the most promising sources for kilometer- 
size laser-interferometric gravitational- wave detectors 
II]. A statistical study based on the stellar evolution 
synthesis (e.g., Refs. suggests that detection rate 

^ 1 — 100 yr~^ may be achieved by advanced detectors 
such as advanced LIGO Q, advanced VIRGO _7], and 
KAGRA (LCGT) [1], which will be in operation in this 
decade. 

One of the key goals after the first detection of gravi- 
tational waves from inspiraling black hole-neutron star 
(BH-NS) binaries and binary neutron stars (NS-NS) 
achieved in the near future will be to extract binary pa- 
rameters such as mass, spin, and radius of each object 
in the binary systems. In particular, the mass and quan- 
tities related to the finite size of the neutron star will 
provide us the essential information for the equation of 
state (EOS) of the neutron-star matter. The mass of two 
neutron stars will be determined with a high accuracy 
< 1%, if the gravitational-wave signals in the inspiral 
stage are detected with the signal-to-noise ratio > 10 9J. 

It is more challenging to determine the parameters 
related to the finite size of neutron stars. Several 
methods to measure such a quantity by observing the 
gravitational-wave signals from NS-NS and BH-NS bi- 
naries have been already proposed fl0l - [T7| . In particu- 
lar, extracting the tidal deformability of neutron stars 
from the gravitational-wave signals from NS-NS inspi- 
rals is one of the convincing ways pTj . For realizing this 
method, one needs to go beyond the point-particle ap- 
proximation to model the gravitational waveform in NS- 
NS inspirals. In other words, one has to derive the grav- 
itational waveform in NS-NS inspirals including tidal ef- 
fects, which influence the dynamics of the binary systems 
in the late inspiral phase (e.g., Ref. (3). When the tidal 
deformability of a neutron star can be measured by the 
gravitational- wave observations, one can constrain on the 



neutron-star matter EOS. Once the neutron-star EOS is 
known, one also can measure the relationship between 
the luminosity distance and the redshift of the binary 
using only the information of the gravitational wave- 
form through the tidal deformation of the inspiralling 
neutron stars [l9j with more refined detectors such as 
Einstein telescope [13]. Therefore, modeling the gravi- 
tational waveform in NS-NS inspirals including tidal ef- 
fects is important not only from astrophysical point of 
view but also from the viewpoint of nuclear physics and 
cosmology. 

For the early stage of NS-NS inspirals (/ < 400 Hz), 
a post Newtonian (PN) gravitational waveform was de- 
rived by Flanagan and Hinderer including the leading- 
order tidal effects [2l| . It shows that the tidal interaction 
affects the evolution of the gravitational- wave phase only 
through a single parameter, namely the tidal deforma- 
bility of a neutron star, up to the leading order. They 
also showed that the tidal deformability of a neutron star 
could be measurable by the advanced gravitational-wave 
detectors by using the gravitational- wave signals for 10 
- 400 Hz, if the tidal deformability of a neutron star is 
sufficiently large or if we observe an event with a high 
signal-to-noise ratio (see also Ref. [22|). 

More recently, focusing on the late stage of NS-NS in- 
spirals (/ > 400 Hz), Damour and his collaborators [23 1 
explored the measurability of the tidal deformability with 
the advanced gravitational- wave detectors. They used 
an effective one body (EOB) formalism for modeling the 
waveform in NS-NS inspirals including tidal effects up to 
2PN order. They concluded that the tidal deformabil- 
ity of a neutron star can be measured by the advanced 
gravitational-wave detectors for the gravitational-wave 
signals of which the signal-to-noise ratio is higher than 16 
for any EOS that satisfies the constraint of the maximum 
mass > l.QTAf© [2^]. The key assumption of their study 
is that the EOB formalism is valid up to the contact point 
of the two neutron stars. 

In the stage just before the merger, nonlinear hydro- 
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dynamics effects play a crucial role for the evolution of 
NS-NS binaries 25| . In addition, higher-PN tidal correc- 
tions may yield a pole-like behavior of the tidal interac- 
tions near the last unstable orbit (26j . 

For better understanding the precise motion and the 
waveform in this late inspiral stage, a numerical-relativity 
(NR) simulation is probably the best approach (see, e.g., 
Refs. (27l - l29j for a review of this field). Recently, long- 
term simulations for NS-NS inspirals were performed by 
three groups (sol - ls^ aiming at the derivation of accu- 
rate gravitational waveforms for the late inspiral stage. 
Baiotti and his collaborators performed a NR simula- 
tion employing a F-law EOS and compared the resulting 
waveforms of the highest resolution simulation with the 
analytic models calculated in the EOB and Taylor T4 
formalisms [s^, HI]. They suggested that the tidal ef- 
fects might be significantly amplified by higher-PN tidal 
corrections even in the early inspiral phase. 

Bernuzzi and his collaborators performed a simulation 
with P-law EOS (P = 2, and the compactness of a neu- 
tron star is 0.14) [3i,|35|. In Ref. "s^, they studied the 
convergence of the numerical results for NS-NS inspirals. 
They concluded that the convergence of the simulation 
is second order up to contact. They also compared the 
resulting extrapolated waveform with that of the Taylor 
T4 formalism for the point-particle approximation and 
for including the tidal corrections. They found that the 
accumulated phase difference is about 1.5 radian at con- 
tact for a particular model of the NS-NS binary. In the 
subsequent paper [ssj . they compared the waveform de- 
rived by the highest resolution simulation with the wave- 
form calculated in the EOB formalism. They found that 
the EOB formalism including tidal corrections up to the 
next-to-next-to leading order is currently the most ro- 
bust way to describe the waveform of NS-NS inspirals. 
In addition, they excluded the huge amplification of the 
tidal corrections suggested in Ref. [soj . 

In this paper, we study NS-NS inspirals by NR sim- 
ulations with three different EOSs and compare the ex- 
trapolated NR waveforms with those calculated in the 
EOB and Taylor T4 formalisms. Here we extrapolate 
NR data with a new extrapolation procedure, the time 
and phase extrapolation. For studying the dependence of 
the tidal effects on the neutron-star matter EOS, we em- 
ploy a piecewise-polytropic EOS of Ref. [36:] , which can 
approximately describe the EOS based on nuclear the- 
oretical calculations and more realistic than F-law EOS 
adopted in Refs. [13, [13, HI] ■ In this paper, (i) we obtain 
the physical gravitational-wave phase by extrapolation; 
(ii) we then compare the extrapolated waveforms with 
those of the analytic models calculated in the Taylor T4 
and EOB formalisms; (iii) we clarify the tidal effects on 
the gravitational-wave phase and show the validity of the 
analytic modeling in the late inspiral phase. 

The paper is organized as follows. In Sec. II, we briefly 
review the analytic modeling of gravitational waves emit- 
ted from a tidally interacting binary system. In Sec. Ill, 
we summarize the formulation and numerical schemes 



employed in our numerical code SACRA, and review the 
EOS employed in this study. In Sec. IV, we describe our 
method of data analysis of the numerical waveforms; the 
radius extrapolation and the resolution extrapolation. In 
Sec. V, we compare the extrapolated gravitational-wave 
phase with those derived in the analytic modeling. Sec- 
tion VI is devoted to a summary. Throughout this paper, 
we adopt the geometrical units of c = G = 1 where c and 
G are the speed of light the gravitational constant re- 
spectively. 



II. TIDAL EFFECTS IN A BINARY SYSTEM 

In this section, we describe analytic models for the 
calculation of gravitational waves emitted from NS-NS 
inspirals in close orbits. We briefly summarize the def- 
inition of the tidal deformability of a neutron star, and 
the PN and EOB descriptions of the tidally-interacting 
dynamics of close NS-NS binaries. 



A. Tidal deformability of a neutron star 

In a close binary system for which the separation be- 
tween two stars is a few times larger than the stellar ra- 
dius, each star is deformed from its hypothetical equilib- 
rium shape in isolation due to the tidal fields. Assuming 
that neutron stars are spherically symmetric in the zeroth 
order, such deformation can be described as the linear re- 
sponses of neutron stars to external tidal fields (37l - [39j . 
as long as the degree of the tidal deformation is small. In 
this linear theory, one assumes that the mass quadrupole 
moment of a star, Qij, is proportional to the external 
quadrupolar tidal fields £ij as, 

Qij — (1) 

where A is the quadrupolar tidal deformability of the star. 
This relation is called the adiabatic approximation for the 
tidal deformation of a star, which is valid only when the 
time scale in the change of the weak tidal field is much 
longer than the dynamical time scale of the star. The 
tidal deformability is related to the quadrupolar tidal 
Love number k2 by 

A-^i?'fc2, (2) 



where R is the radius of the (spherical) star in isolation. 
For a given EOS and a central density, one can calculate 
the quantities mass, R, k2, and A of neutron stars by 
solving the Tolman-Oppenheimer-Volkoff equations and 
the metric perturbation equations [13, HI] • 



B. The post- Newtonian description for the motion 
of a tidally interacting binary 

The motion of tidally interacting NS-NS binaries in 
close orbits is affected by the stellar internal structure. 
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As long as the degree of the tidal interaction is small, the 
correction of this effect can be described only through the 
tidal deformability A [i^. The evolution of the orbital 
angular velocity in the inspiral of a tidally interacting 
binary is described by 



NS-NS inspirals only at 5PN order, its coefficient is of 
order 10'* for typical neutron stars of radius 10 - 15 km, 
/c2 ^ 0.1, and C ^ 0.14 - 0.20, and thus, it plays an 
important role in the late inspiral stage. 



dx 
H 



F{x,Ma,Mb,Xa,\b), 



(3) 



where x = (Muj)'^/^ with uj being the angular velocity 
of the binary and M being the total ADM mass at the 
infinite separation. The subscript (A or B) refers to each 
component of the binary. The function F can be decom- 
posed into the following two parts, 

F = Fpp(x, Ma, Mb) + Fud^i{x, Ma, Mb, Xa, \b), (4) 

where _Fpp is the contribution of the point-particle part, 
and -Ftidai is the contribution associated with the tidal 
interactions. In this work, we adopt the Taylor T4 ap- 
proxiniant for the i^pp [4l| . 
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where — MaMb/M"^ is the symmetric mass ratio, and 
7 = 0.577126... is an Euler's constant. As shown by Boyle 
and his collaborators |4l|, the evolution of the angular 
velocity of the Taylor T4 approximant agrees well with 
those of the NR simulations for the inspiral of equal- 
mass non-spinning binary black holes up to Mu 0.1. 
We adopt the tidal part which is derived by Vines and 
his collaborators yol as follows, 
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In particular, for the case of equal-mass binaries, Ftidai 
is given by 



52 ^ 5 A 5203 

Ftidai = x^'^k2C^'' 1 H X 

5M ^ I 4368 



(7) 



C. Effective one body formalism for the motion of 
a tidally interacting binary 

The EOB formalism maps the dynamics of two point 
particles to the Hamiltonian dynamics of an effective par- 
ticle moving in an effective external potential |43l - |45| . 
Because the EOB formalism goes beyond the adiabatic 
approximation for binary inspirals, it is suitable for de- 
scribing the late inspiral stage of a binary, for which the 
adiabatic approximation is not very accurate. In this 
work, we employ the resummed EOB description which is 
largely the same as that in Ref. [i^ for the point-particle 
part and as that in Ref. 23] for the tidal part. 

The EOB effective metric is defined by 



dsl, = -A{r)dt' + ^dr' + 
A{r) 



+ sin^ 0d<f) , (8) 



where (t, r, </)) are the dimensionless coordinates and their 
canonical momenta are (j)r,Pcj,)- We replace the radial 
canonical momentum pr with the canonical momentum 
Pr , where the tortise-like radial coordinate is defined 

by' 



dr 



A{r) ■ 



(9) 



Then the binary dynamics can be described by the EOB 
Hamiltonian 



i?rcal(r,Pr.,P0) = ^yl + 2l/(ifcff-l) (10) 

where the effective Hamiltonian is defined by. 



Hcff — 



Pi 



A(r) 1 + ^ + 2(4-3^.)^.^ . (11) 



The metric component A{r) is decomposed into two parts 



A{r) = App(r)+Atidai(0, 



(12) 



where App{r) is the point-particle part and ^ItidaK?") is 
associated with the tidal effects. The point-particle part 
up to the 5PN order is given by 



^PP('^) = P5 



1 — 2m -|- 2i>u^ + a4iyu'^ 



, (13) 



where C = Ma/ Ra {= M-q/ Rb) is the compactness of 
the neutron star. Although the tidal interaction affects 



where 04 = 94/3 — 4l7r^/32, u = 1/r, and denotes a 
(1, 5) Fade approximant. In the definition of A{r), there 
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are two analytically undetermined parameters (05, ag), 
which correspond to the 4PN and 5PN corrections. Here, 
we adopt the values (05,05) = ((—5.828 — 143. + 
477z/2)jy, 184iy) following Ref. 

The tidal part of A{r) is given by 



A 



kidal 



(14) 



l>2 



where Ai includes the PN tidal contributions for each 
multipole, and ki is its coefficient. This coefficient is 
related to the tidal Love number ki and the compactness 
of two stars by 



21 



Kl = 2 



kf 



M21+1 ^2/4 



(A^B) 



(15) 



In this work, we include only the tidal-interaction part 
of the lowest multipole 1 = 2. Up to the next-to-leading 
order, Ai=2 is given by [42j 



(16) 



where a^*^ = 5Ma/2M. The tidal-interaction term up 
to the 2PN corrections is currently known, and the coeffi- 
cient a2 is larger than ai '26']. In addition, an analysis in 
the test- mass limit (Ma ^ Mb) suggests that the tidal 
part ^1^2 the pole-like behavior near the last unsta- 
ble orbit located at 3Mb (the light ring orbit). Thus, 
it is reasonable to expect that the higher-PN corrections 
would amplify the tidal effects. In this work, we employ 
the resummed version of the tidal metric including up to 
the next-to-next-to-leading order given by [26] 



= 1 



1 - ^lru ' 



where 



and 



a. 



(A) 



337Afl/28Af2 + Ma/SM + 3, 



rhRiv,K2) = 3 



4 2 2n 

33 ^ ^'«2 + 0(V , K2l^, K2) 



(17) 



(18) 



, (19) 



is the dimensionless radius of the light ring orbit. 

For the calculation of the binary orbit, we solve the 
BOB Hamihon equations psj-ISsj 



dr 
'dt 



A{r) dH, 



real 



dpr, 

dt 



d4 _ giJrcal 
dt ~ 

A{t) 



real 



dp,t> ir 
-dt=^''' 



(20) 
(21) 
(22) 
(23) 



where D^{r) is a (0,3) Pade approximant of D{r) 
and Jv, is the radiation-reaction force given by 



1 



1=2 m=l 



M 



(24) 



Here, u) = d(j)/dt, him denote the multipolar waveforms, 
and R is the radius of extracting gravitational waves. 
The waveforms are described by 



Hm 



'hm 



- tidal, A 



1 tidal,; 



(25) 



where hf^ denotes the inspiral and plunge waveform for 



J 7 tidal. A . 

and /i,^ IS 



a binary black hole of mass Ma and Mb 
the contribution due to the tidal deformation of star A 
In this work, /i?.^ is given by Eqs. (13)-(22) of Ref. [i 



and ^''^ is basically given by Eqs. (A14)-(A17) of 
Ref. [23|. Because the 2PN term associated with the 
tidal effect in the waveform is currently unknown, we 



include the contributions of the tidal effect to h 



tidal, A 
Im 



up 



to the IPN term even for the case that the next-to-next- 
to-leading order term is taken into account in the radial 
potential. 



III. NUMERICAL RELATIVITY SIMULATION 

In this section, we briefly describe the formulation and 
the numerical schemes of our NR simulation employed in 
this work. 



A. Evolution and Initial Condition 

We follow the inspiral and merger of NS-NS bina- 
ries using our NR code, called SACRA, for which the de- 
tails are described in Ref. [47]. SACRA employs a mov- 
ing puncture version of the Baumgarte-Shapiro-Shibata- 
Nakamura formalism [isl - fsoj to solve Einstein's equations 
imposing the equatorial symmetry (and 7r-symmetry for 
the equal-mass cases). In the numerical simulations, a 
fourth-order finite differencing scheme in space and time 
is used implementing an adaptive mesh refinement algo- 
rithm. At refinement boundaries, a second-order inter- 
polation scheme is partly adopted. The advection terms 
are evaluated by a fourth-order non-centered finite dif- 
ferencing [5l|. A fourth-order Runge-Kutta method is 
employed for the time evolution. For the hydrodynam- 
ics, we employ a high-resolution central scheme based 
on Kurganov and Tadmor scheme [s^ ] with a third-order 
piecewise parabolic interpolation and with a steep min- 
mod limiter. 

In this work, we prepare seven refinement levels both 
for efficiently resolving two neutron stars by the finest- 
resolution domains and for extracting gravitational waves 
in a local wave zone. More precisely, two sets of four 
finer domains comoving with each neutron star cover 



5 



TABLE I: Key parameters for the initial models adopted in the present numerical simulation. Mo is the sum of the ADM 
masses of two neutron stars in isolation; u is the symmetric mass ratio; Mq^^ and J^^^ are the ADM mass and angular 
momentum of the system, respectively; Af, is the baryon rest mass; ujq is the angular velocity. We also show the setup of the 
grid structure of our AMR algorithm. Aa; — ha = L/{2^N) {N = 60) is the grid spacing for the highest-resolution domain with 
L being the location of the outer boundaries for each axis. A'^ specifies the grid size of the simulation with maximum N = 60. 
K2 is the parameter related to the tidal deformability. tJcontact is the orbital angular velocity at contact of the two neutron stars 
(see Sec. V. B). Here we use the unit Mq = 1 



Model 


M, 


u 




tAUM 


M. 


Mojo 


Ax/M 


Ii2 


A'/oJcontact 


N 


APR4 


2.7 


0.25 


2.68 


7.65 


3.00 


0.019 


0.0438 


62.3 


0.151 


(42,48,54,60) 


H4 


2.7 


0.25 


2.68 


7.66 


2.94 


0.019 


0.0560 


215 


0.112 


(42,48,54,60) 


MSI 


2.7 


0.25 


2.68 


7.67 


2.92 


0.019 


0.0595 


332 


0.103 


(42,48,54,60) 


APR4-1215 


2.7 


0.24691358 


2.68 


7.28 


3.01 


0.0221 


0.0438 


65.8 


0.151 


(40,50,60) 


H4-1215 


2.7 


0.24691358 


2.68 


7.56 


2.94 


0.019 


0.0560 


207 


0.114 


(40,50,60) 



TABLE IL Parameters of the piecewise-polytropic EOS, the compactness, and the tidal Love number of the neutron star of 
mass 1.35M0. 



EOS 


log Pi (dyne/cm^) Pi Pa P3 C(1.35) fc2(1.35) 


APR4 
H4 
MSI 


34.269 2.830 3.445 3.348 0.179 0.091 
34.669 2.909 2.246 2.144 0.146 0.115 
34.858 3.224 3.033 1.325 0.138 0.133 



the region of their vicinity. The other three coarser 
domains cover both neutron stars by a wider domain 
with their origins fixed approximately at the center of 
the mass of the binary. Each refinement domain con- 
sists of the uniform, vertex-centered Cartesian grids with 
{2N + 1,2N +1,N + 1) grid points for {x, y, z) with the 
equatorial plane symmetry at z = imposed. The half 
of the edge length of the largest domain (i.e., the dis- 
tance from the origin to outer boundaries along each axis) 
is denoted by L which is chosen to be larger than Ao, 
where Ao = tt/wq is the initial wave-length of gravita- 
tional waves and ljq is the initial orbital angular velocity. 
The grid spacing for each domain is hi — L/{2^N), where 
I — Q — Q. In this work, we choose N = 60, 54, 48, and 
42 for the resolution study. With the highest grid res- 
olution, the semimajor diameter of each neutron star is 
covered by about 100 grid points. 



We prepare NS-NS binaries in quasiequilibrium states 
for the initial condition of numerical simulations by using 
a spectral- method library, LORENE [5I] . To track more 
than 8 orbits, the orbital angular velocity of the initial 
configuration is chosen to be Mwq = 0.019 (/ = 400Hz 
for M = 2.7 Mq), where Mq is the solar mass. The 
neutron stars are assumed to have an irrotational veloc- 
ity field, which is believe to be an astrophysically realis- 
tic configuration [HJ, El] • The parameters for the initial 
models are listed in Table HI 



B. Equation of State 

In this work, we employ a parameterized piecewise- 
polytropic EOS proposed by Read and her collabora- 
tors [36,] . This EOS is written in terms of four segments 
of polytropes 

P ^K,p^' (26) 
( for Pi< p< Pi+i, < i < 3), 

where p is the rest-mass density, P is the pressure, 
Ki is a poly tropic constant, and Ti is an adiabatic in- 
dex. At each boundary of the piecewise polytropes, 
p = Pi, the pressure is required to be continuous, i.e., 
Kip^^ = Ki+ip^"''^. Read and her collaborators de- 
termine these parameters in the following manner (36|. 
First, they fix the EOS of the crust as Tq = 1.357 and 
Kq ~ 3.594 X 10^"^ in cgs units. Then they determine 
P2 = 1.85pnuci and ps = 3.70(Onuci where pnud = 2.7x 10^^ 
g/cm"^ is the nuclear saturation density. With this prepa- 
ration, they choose the following four parameters as a 
set of free parameters: {Pi, Fi, Fa}. Here Pi is the 
pressure at p = p2, and from this, Ki and Ki are deter- 
mined by Ki = P\l ^^'^ Ki^i = Kip^' '"•+\ There- 
fore the EOS is specified by choosing the four parameters 
{Pi, Fi, F2, F3}. In this work, we choose the three sets of 
piecewise-polytropic EOS as listed in Table |lll 

We describe the low-density part of the EOS only with 
a single polytrope, because the elastic property of the 
crust yields a very small correction to the tidal num- 
ber [58|. Thus, our approximate treatment for the low- 
density part is acceptable. 
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C. Extraction of Gravitational waves 

Gravitational waves are extracted by calculating the 
complex Weyl scalar ^'4 from which gravitational 

waveforms are determined by 

h+{t) -ih^{t) = - \im [ dt' [ dt"^i{t'\r). (27) 

Here we omit arguments and (j). In the spherical coor- 
dinate {r,6,(f>), ^'4 can be expanded in the form 

^^{t, r, 0,<P)=Y. *4'"(^, r)^2Yi„,{0, 0), (28) 

Im 

where -2Yim are spin- weighted spherical harmonics of 
weight —2 and ^'4" are expansion coefficients defined 
by this equation. In this work, we focus only on the 
{I, |m|) = (2, 2) mode. The gravitational- wave phase (/)nr 
is defined by 

*f (i,r) = A22(t,r)e*^~«(*''^), (29) 

where A22 denotes the amplitude and it is real. We eval- 
uate ^'4 at a finite spherical-coordinate radius t/Mq — 
200, 240, 300, and 400. To compare the waveforms ex- 
tracted at different radii, we use the retarded time defined 
by 

irot = t-r^. (30) 
Here, is the tortoise coordinate defined by 

=rA + 2A/ln(^-l), (31) 

where ta = y/A/An with A being the proper surface area 
of the extraction sphere. 

IV. DATA ANALYSIS OF 
NUMERICAL-RELATIVITY SIMULATIONS 

The extrapolation with respect to the extraction radius 
and grid resolution plays a key role to obtain the physical 
waveforms in NS-NS inspirals from the results of NR sim- 
ulations. Here, we focus in particular on deriving the ac- 
curate gravitational- wave phase by extrapolation because 
it carries the most important information in the matched 
filtering for data analysis. Thus, the goal of this section is 
to construct the extrapolated gravitational- wave phase. 

A. Extrapolation to infinity 

Because gravitational waves are extracted at finite 
radii, the extracted waveform does not fully agree with 
the waveform at infinity. For obtaining the hypothetical 
gravitational- wave phase at infinity, we first need to esti- 
mate the error due to the finite-radii extraction and then 



to extrapolate the gravitational-wave phase to infinity. 
For this purpose, we assume that t he g ravitational-wave 
phase is described by a polynomial [41|, 

0(t,,„r) = </-(o)(i,,O+E^^^- (32) 

i=l 

Here, 0^°^(irct) is considered to be the gravitational- wave 
phase extrapolated at r — 00, and (s -I- 1) is the num- 
ber of extraction radii used. We determine it by ex- 
trapolating the gravitational-wave phases extracted at 
t/Mq = 200, 240, 300, and 400. 

Figure [1] shows the differences among the hypotheti- 
cally extrapolated gravitational-wave phases at infinity 
obtained by different numbers of extrapolation radius, s. 
The curves labeled by 1 - 3(= s — 1) in Fig. [1] denote the 
differences between a radius-extrapolated gravitational- 
wave phase with the (s — l)-th order polynomial and 
that with the s-th order one. Here, we set the phase 
difference to be zero at the merger. Note that the dif- 
ference between the (s — l)-th and s-th order extrapo- 
lated gravitational-wave phases accumulates mainly just 
after the contact of the two stars and its value is less 
than about 0.5 radian. In this paper, we do not pay 
attention to the gravitational-wave phase after contact 
of the two stars. As shown in the next subsection, 
furthermore, this accumulated gravitational-wave phase 
difference of the radius extrapolation is much smaller 
than the magnitude of other error sources. Therefore, 
in this paper, we neglect the phase error associated 
with the finite-radius extraction, and use the resolution- 
extrapolated gravitational- wave phase with the extracted 
radius t/Mq = 400 as the extrapolated wave phase. We 
note however that in a high-resolution study, the error 
associated with the finite-radius extraction could be com- 
parable to or smaller than the error due to the finite grid 
resolution. 



B. Extrapolation to infinite resolution 

A NR waveform is derived as a solution of discretized 
Einstein's equations and hydrodynamics equations. As- 
sociated with the finite differencing, a truncation error 
is yielded. For the purpose of comparing the NR wave- 
form with the FN and EOB waveforms, we have to obtain 
a hypothetically physical waveform by extrapolating the 
NR waveform obtained in the finite grid resolution for 
the limit of the infinite grid resolution. The order of the 
convergence of our NR simulation with respect to the 
grid spacing is found to be about 1.8. Because the order 
of the convergence has an uncertainty, we conservatively 
assume that the order of the convergence is 1.8 ± 0.2 in 
the following. 

For understanding the convergence property of the nu- 
merical results, we check the convergence for the evolu- 
tion of the gravitational-wave frequency performing sim- 
ulations with four different grid resolutions (see Fig. [5] 
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FIG. 1: The phase errors for APR4 (left) and H4 (right). The curve 1 denotes the difference between the gravitational- 
wave phase extracted at r/M© — 400 and an extrapolated gravitational-wave phase obtained using the data at t/Mq = 300 
and 400. The curve 2 denotes the difference between the extrapolated gravitational-wave phase obtained using the data at 
r/M© = 300 and 400 and that obtained using the data at r/M© — 240, 300, and 400. The curve 3 denotes the difference 
between the extrapolated gravitational-wave phase obtained using the data at v/Mq = 240, 300, and 400 and that obtained 
using the data at r/M© = 200, 240, 300, and 400. 
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FIG. 2: The evolution of gravitational-wave angular velocity for four different resolutions for APR4 (left panel) and H4 (right 
panel). The curves with marks A, B, C, and D denote the results for the simulations with = (42, 48, 54, 60), respectively. 



for gravitational- wave angular velocity). Here, we re- 
fer to the data for N = (42,48,54,60) as (A, B, C, D). 
Although the gravitational-wave frequencies agree with 
each other at t-^ct ~ 0, their subsequent evolution dis- 
agrees with each other. Obviously, the gravitational- 
wave frequency in the lower grid resolution evolves more 
rapidly than that in the higher grid resolution. The rapid 
evolution for the lower-resolution simulations may be as- 
cribed to larger numerical dissipation of the angular mo- 
mentum. 



However, one can find the similarity among four curves 
of the gravitational- wave frequency for the different grid 
resolutions and this enables us to obtain an extrapolated 
waveform. To show this similarity, we normalize the time 
variable by using the time at the onset of the merger, 



which is defined by tm{^x) — t\^^^,^^{lS.x): 
t i{t) = 



A(Ax) ■ 



(33) 



Here, A(Aa;) — fm(0)/im(Aa;), Aa; is the grid spacing of 
the simulation (oc N~^), and im(0) is the extrapolated 
merger time. For obtaining the extrapolated merger 
time, we assume that the merger time as a function of 
grid resolutions is described by a binomial 



U(Aa;) =t„(0)+X(Aa;)", 



(34) 



where tm(0), ?^, and K are constants which are deter- 
mined by the least-square fitting method. We find that 
the best fitted value of n is « 1.8 with a dispersion 
~ 0.2, and thus, we set the order of the convergence 
n is to be 1.8 ± 0.2. The resulting fitted values for 
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FIG. 3: Convergence of the rescaled gravitational-wave angular velocity. Upper panels: the evolution curves of the rescaled 
gravitational- wave angular velocity for APR4 (left panel) and H4 (right panel). Lower panels: the difference in the rescaled 
gravitational-wave frequency between different resolutions for APR4 (left panel) and H4 (right panel). Here we choose n = 1.8. 



n = (1.6,1.8,2.0) are t™(0) = (2521M, 2452M, 2397M) 
for APR4, i„(0) = (2391M, 2274M, 2238M) for H4, and 
t„, = (2145M,2145M,2091M) for MSI. 

We define the rescaled gravitational- wave frequency by 

Cj{t,Ax) =uj{i{t),Ax). (35) 

The evolution curves of the rescaled gravitational-wave 
frequency u for the different grid resolutions agree with 
each other, as shown in the upper two panels of Fig. [3] 
(here, we set n = 1.8). The differences among the 
rescaled gravitational-wave frequencies of the different 
grid resolutions are within ~ 10% for APR4 and ^ 5% for 
H4, as shown in the lower panels of Fig. [3] Therefore, we 
conclude that the numerical results for the gravitational- 
wave frequency (angular velocity) is written in the form 

w(t. Ax) = w(t,0) Ax), (36) 

where Wi.(i, Ax) is a function that depends on the grid 
resolution, which should have been eliminated system- 
atically by extrapolation. However, because the value of 
Uj-it, Ax) randomly fluctuates, it cannot be eliminated by 



extrapolation at each given moment. Thus, the extrapo- 
lated gravitational-wave frequency can be only approxi- 
mately obtained by simply rescaling the time variable of 
the simulations. 

We proceed to extrapolate the rescaled gravitational- 
wave phase, which will be compared with the 
gravitational-wave phases derived in the PN and EOB 
approaches in the next section. We define the rescaled 
gravitational-wave phase as 

4>{t,Ax) EE [ Cj{t' , Ax)dt' . (37) 
Jo 

Substituting Eq. (|55|) into this equation yields 

^(i. Ax) = A(Ax)(?!)(i, Ax), (38) 
and combining Eq. ([55]) and Eq. ([57)) yields 

(t){t,Ax)^4){t,0)+ [ t(t',Ax)dt', (39) 

where 4)(t, 0) is the extrapolated gravitational-wave 
phase. Assuming the second term in the right-hand side 
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FIG. 4: The extrapolated gravitational-wave phase for the 
three EOSs. Here, we set the order of the convergence to be 
1.8. 



of this equation can be described as 4>i{t)Ax'^, we obtain 

A{Ax)(j){i, Ax) = 4>{t, 0) + (40) 

Here 4>{t,0) and are functions which are deter- 

mined by the least-square fitting of the numerical data 
(/)NR(i, Ax). We again set the order of the convergence n 
to be 1.8 ± 0.2. Figure m shows the evolution curves of 
the extrapolated gravitational-wave phase for the differ- 
ent EOS. 

The differences in the rescaled gravitational-wave 
phase 0nr(^, Ax) obtained for runs with the differ- 
ent grid resolutions and the extrapolated gravitational- 
wave phase (l){t,0) for n — 1.8 are shown in Fig. m 
The accumulated difference between the extrapolated 
gravitational- wave phase and that of the highest grid res- 
olution at the onset of the merger is '-^ 8 radian for APR4 
and 3 radian for H4. These values imply that the sim- 
ulation for the neutron star which has the larger radius 
is relatively well convergent. These phase differences are 
much larger than the possible error in the gravitational- 
wave phase associated with the finite-radius extraction, 
< 0.5 radian. Thus, we conclude that the resolution ex- 
trapolation of the gravitational- wave phase is essential to 
construct the extrapolated physical waveform for our cur- 
rent NR simulations. The lower panels in Fig. [5] show the 
difference between the gravitational-wave phase of the 
resolution D, (^nr(^, Axu), and the one which is obtained 
by extrapolating the gravitational- wave phase of the res- 
olutions A, B, and C to the resolution D, fp{t, Axu)- The 
difference A0 is within 0.15 radian up to the merger. 
Here, this error is caused mainly by a modulation associ- 
ated with an unphysical orbital eccentricity. The upper 
panels of Fig.[6]show the evolution curves of the extrapo- 
lated gravitational- wave phase for the different choice of 
n. 

For the brief check of the validity of the extrapolation, 



we calculate the dispersion a{t) of the extrapolation as 

1 / \2 
'^W' = — ^(0(t,Ax,)-^NR(i,AxOj , (41) 

where (/>NR(i, Axi) is the rescaled gravitational- wave 
phase of the numerical data, i = 1 TVr denote 
the numerical run with the different grid resolutions, 
e.g. N = (42,48,54,60), and TVr is the total number 
of them, e.g. TVr = 4. As shown later, the value of the 
dispersion is in the range 0.01 radian (in the early part) 
to 0.4 radian (just before the merger). We regard this 
dispersion as an error due to the resolution extrapola- 
tion. Hereafter we use only the extrapolated quantities 
and omit the tilde of them. 



V. COMPARISON OF THE NR WAVE PHASE 
WITH THE PN AND EOB WAVE PHASES 

A. Matching Procedure 

In the analysis, we first match an early part of the 
gravitational-wave phase obtained by extrapolating NR 
results , 0NR(i), with that of the analytic ones, (/)pN(i) 
(the wave phase derived in the PN calculation) and 
0EOB {t) (the wave phase derived in the EOB calculation) 
, by minimizing the following quantity, 

f\A<f>it))^ dt 




(0NR(i) - '/'PN/EOB(i - is) - 0s) dt, (42) 



where (ps and tg are the fitting parameters. The initial 
and final time of the integral, ti and t2, are chosen as fol- 
lows. As shown in the previous section, the curve of the 
extrapolated gravitational- wave phase obtained from the 
results of NR simulations has a small modulation due to 
the orbital eccentricity. Thus, (ii,i2) are chosen to cover 
a range between two adjacent local maxima or local min- 
ima of the curve in an early part of the gravitational- 
wave phase. This choice of the matching region al- 
lows us to match the gravitational-wave phase numeri- 
cally obtained with the PN/EOB phases smoothly using 
the least-square fitting. Here, we choose the two adja- 
cent local minima as the matching boundaries as follows: 
{ti,t2) = (583M, 1030M) for APR4 (which corresponds 
to Moj e [0.042 : 0.048]), (ti,t2) = (550M, lOOOM) for 
H4 (which corresponds to Muj e [0.042 : 0.046]), and 
(ti,t2) = (550M,960Af) for MSI. 

B. Comparison 

We compare (/)nr with ippjq and (^eob which are calcu- 
lated in several levels of the approximation. The middle 
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FIG. 5: Convergence of the gravitational- wave phase for APR4 (left) and H4 (right). In the upper panels, the curves show the 
phase differences between the extrapolated gravitational-wave phase and the rescaled numerical results. Here, we set the order 
of the convergence n to be 1.8 for deriving the extrapolated gravitational-wave phase. In the lower panels, the curves show 
the difference between the rescaled gravitational-wave phase in the resolution D and the extrapolated gravitational-wave phase 
which is obtained by extrapolating the rescaled gravitational-wave phases of the resolution A, B, and C to the resolution D. 



panels of Fig. [5] show the difference between (/)nr and 
0EOB in which tidal effects up to the next-to-leading or- 
der are included. The error bars show the dispersion as- 
sociated with the least-square fitting in the extrapolation 
procedure, which is defined by Eq. (41). The three curves 
show the gravitational- wave phase difference for the three 
assumed orders of the convergence n = (1.6,1.8,2.0). 
The phase differences are modulated with an amplitude 
about 0.4 radian due to the unphysical orbital eccentric- 
ity. We regard this modulation as a systematic error, 
which is denoted by the horizontal dashed lines. For 
comparison, we also plot the gravitational-wave angular 
velocity evolution in the bottom of Fig. [B] Figure [7HS1 
shows the difference between (^nr and (?!ipN/EOB which 
are calculated in several levels of the approximation. 

In the early stage of the inspiral, i.e., before the time 
tret/M < 1500 (Mw < 0.054) for APR4, Uet/M < 
1300 [Muj < 0.05) for H4, and Uct/M < 1200 (Mw < 
0.049) for MSI, we find that ^nr is consistent with (pp-^ 
and 4>EOB for any levels of the approximation. In this 
stage, the point-particle approximation works well and 
the FN and EOB approaches appear to describe NS-NS 



inspirals well even if tidal effects are not taken into ac- 
count. It is difficult to verify a clear signature of the tidal 
effects in this early stage due to the modulation of the 
numerical data and the weakness of the tidal effects. 

After the early stage of the inspiral, the binary sys- 
tem proceeds to a tidally-dominated inspiral phase where 
the tidal interaction becomes strong. In this late inspi- 
ral stage, one can see that the difference between (/)nr 
and (/ipN/EOB gradually increases. 0nr, is larger than 
0PN/EOB for a given moment, implying binaries in NR 
simulations evolve faster than those in the FN and EOB 
calculations. 

After the time t.^t/M - 2300 {Moj - 0.09) for APR4, 
trct/M - 2100 {Muj - 0.08) for H4, and Uet/M ~ 
f 950 {Muj - 0.075) for MSI, the difference between (/.nr 
and (fiEOB rapidly increases with increasing time. This 
corresponds to the transition from the inspiral to the 
plunge. The orbital angular velocity when two neutron 
stars come into contact, Mwcontact, is approximately de- 
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At the moment of this contact, the accumulated 
gravitational-wave phase difference between (/)nr and 
0EOB including the tidal effects up to the next-to-leading 
order is ^ 3.3 radian (2.5%) for APR4, ^ 1.9 radian 
(1.7%) for H4, and - 2.1 radian (2.0%) for MSI. If 
the correction is up to the next-to-next-to-leading order, 
the phase difference is ~ 2.6 radian (2.0%) for APR4, 
- 1.4 radian (1.3%) for H4, and - 1.1 radian (1.1%) for 
MSI. We find that the EOB approach including the tidal 
effects up to the next-to-next-to-leading order yields cur- 
rently the best model for the late stage of NS-NS inspi- 
rals. However, the tidal effects are still underestimated 
for the final inspiral orbit even for the best model of 
the analytic approaches. Thus, for constructing better 
waveforms in NS-NS inspirals, this rapid evolution of the 
gravitational-wave phase should be taken into account. 

We compare our results of MSI, of which the com- 
pactness is C = 0.138, with the results of Ref. 3J|, for 
which a neutron-star with C — 0.14 is employed. Both 
results show that the difference between the extrapolated 
gravitational-wave phase of NR and that of T4 without 
tidal effects is about 3 radian at contact [Mljj ^ 0.1), in 
the case that the alignment is performed after the first 
orbit [131 • Therefore our results are consistent with the 
results of Ref. fs^. 

Figure [TUl and Fig.fTTjshow the snapshots of the density 
contour of NS-NS inspirals for APR4 and H4. Here we 
focus on Fig. [TT] as an example. The upper left panel 
of Fig. [TT] plots the density profile of the binary at an 
early inspiral stage. At this time, the neutron stars have 
a spherical shape. The ellipticity of the neutron star, 
which is defined by the ratio of the semi-major axis ai 
to the semi-minor axis 02 of the star on the equatorial 
plane, is approximately unity. 

The upper right panel of Fig.[TTl at — 0.056, plots 
the configuration of the binary at which the tidal effects 
seem to be small but cannot be neglected. For this plot, 
the ellipticity of the star is ~ 1.05. After this time, the 
binary system proceeds to the tidally-dominated inspi- 
ral phase. In the lower left panel of Fig. [TTJ two neu- 
tron stars are obviously deformed due to the strong tidal 
fields; the ellipticity of the neutron star is ~ 1.17. The 
snapshot around the plunge is shown in the lower right 
panel of Fig. [TTJ Soon after the onset of the plunge, the 
two neutron stars contact and the ellipticity is ~ 1.23. 
In addition, one can see the appearance of the dynamical 
tidal lag. It appears even in the absence of viscous dis- 
sipation, because the shape of the star cannot follow the 
rapid change of the tidal potential (see, e.g., Ref. [soj). 
Therefore, the adiabatic approximation for the tidally in- 
duced quadrepole moment determined by Eq. (1) breaks 
down after the appearance of the dynamical tidal lag. In 
such a case, the tidal deformability should be evaluated 



with the formalism beyond the adiabatic approximation, 
which is formulated in Ref. [60| . 

We also note that the ellipticity of the neutron star 
rapidly increases after the stage at Muj = 0.056. This 
evolution of the ellipticity is consistent with the semian- 
alytic, Newtonian results [25|. Note that the value of the 
ellipticity defined here depends on the coordinate system. 
We here assumed that the coordinate distortion is small 
because of our choice of the spatial gauge condition. 

Finally, we show the difference between the extrap- 
olated gravitational-wave phase and the gravitational- 
wave phase derived in the FN and EOB approaches 
for unequal- mass systems APR4-1215 and H4-1215: 
see Fig. [T^] . Here, we choose the matching region 
(ii,t2) = (520M,1020M) for APR4-1215 and (^1,^2) = 
(600M, 1080M) for H4-1215. This figure shows that the 
feature of the curves for the unequal- mass system is qual- 
itatively the same as those of the equal-mass system, and 
the magnitude of the phase difference between (/)nr and 
'/'PN /EOB is also as large as that for the equal-mass case. 



VI. CONCLUSION 

We have explored the property of gravitational waves 
emitted in the late stage of NS-NS inspirals. To derive the 
physical waveform from numerical data obtained by NR 
simulations, we carefully performed the radius and reso- 
lution extrapolation of the waveforms. Then, we found 
the resolution extrapolation is crucial in our present 
study. Specifically, the accumulated difference between 
the gravitational-wave phase in the highest grid resolu- 
tion run and the resolution-extrapolated gravitational- 
wave phase is ^ 8 radian for APR4 and ^ 3 radian for H4. 
These values imply that the simulation for the NS-NS 
inspiral with more compact neutron stars has worse con- 
vergence. For the simulation with more compact neutron 
stars , one needs to perform a higher resolution simula- 
tion to derive an accurate waveform. Therefore a sophis- 
ticated procedure for the extrapolation is needed to de- 
rive an accurate waveform for compact neutron stars. We 
found that the time-rescaling is a robust prescription for 
deriving the resolution-extrapolated gravitational-wave 
phase </)NR- 

We have compared (/)nr with (/)pN/EOBi which are de- 
rived from the Taylor T4 approximant of the PN for- 
malism and the EOB formalism. Both of the analytic 
approaches are capable of including the tidal effects. We 
found that 0nr is consistent with (/>pn/eob in the early 
part of the inspiral. On the other hand, in the very 
late part of the inspiral, 0nr evolves more rapidly than 
(/)pN/EOB- The EOB approach including the tidal cor- 
rections up to the next-to-next-to-leading order is cur- 
rently the best approach for describing the late stage of 
NS-NS inspirals. However, the estimated accumulated 
difference between (/)nr and (/!)eob is ~ 2.6 radian for 
APR4, - 1.9 radian for H4, and - 1.1 radian for MSI 
at the moment of contact of the two neutron stars. We 
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FIG. 6: Gravitational- wave phase (top panels), diflFerence between 0nr and 0eob (middle panels), and gravitational- wave 
angular velocity (bottom panels) for APR4 (left panels) and H4 (right panels). EOB(NLO) denotes the effective-one-body 
formalism including the tidal effects up to the next-to leading order. We show the curves for the three orders of the convergence, 
n = (1.6,1.8,2.0). The horizontal dashed lines in the middle panels denote the uncertainty due to the modulation of the 
numerical data. The vertical dashed lines denote the time at contact of the two stars for the case of the order of the convergence 
is 1.8. 



conclude that the tidal effects are still underestimated 
in the EOB approach including the tidal corrections up 
to the next-to-next-to-leading order. We also found that 
this result is independent of the EOS and mass ratio. 

Here we make a comparison of our results with the ear- 
lier results. We find the absence of the large amplification 
of the tidal effects in the early inspiral phase suggested in 
Ref. [30] , if the resolution extrapolation is taken into ac- 
count. In the late inspiral phase, the amplification of the 
tidal effects is observed. This agrees with the earlier re- 
sults HO, [33. In particular, for the MSI (C ~ 0.14), the 
value of the phase difference between the extrapolated 
gravitational-wave phase and that of Taylor T4 without 



tidal effects is consistent with the results of Ref. [3^ . 
which employed the neutron star with C = 0.14. 

For extracting the tidal deformability of a neutron star 
efficiently and faithfully from a signal of gravitational 
waves, one has to prepare a theoretical template of grav- 
itational waves which should be accurate enough up to 
the onset of the merger. Our present study suggests that 
the EOB approach including up to the next-to-next-to- 
leading order tidal correction yields currently the best 
result. However, for the final orbit, there is still room for 
the improvement. In the current prescription for incor- 
porating the tidal effects, the adiabatic approximation is 
assumed for the tidal deformation. For the very close 
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orbits, however, this approximation breaks down; for ex- 
ample, the presence of the dynamical tidal lag (which is 
seen in Figs. ITOl and fTTj) cannot be reproduced. If a more 
sophisticated formalism in which such effects 

4 I ^ ^ ^ ^ T-T, 

NR-T4(pp) i, 

NR-EOB(pp) 

NR-T4(NL0) - ; 

3 - NR-EOB(NLO) I - 

NR-EOB(NNLO) - - - j 




FIG. 7: Phase difference between (^nr andf/ipN/EOB for 
APR4. In this figure, we set the order of the convergence 
to be 1.8 for obtaining the extrapolated gravitational-wave 
phase. T4(pp) and EOB(pp) denote the Taylor T4 approxi- 
mant and the effective-one-body formalism without the tidal 
corrections. Here, we only display the error bars of the curve 
NR-EOB(NLO). However, other curves also have the same 
amount of the error. 



are taken in account can be developed, the accuracy 
of the analytic modeling could be improved (see e.g., 
Ref 

There is also the issue on the side of numerical rela- 
tivity. For constructing an analytic model of the gravita- 
tional waveform in NS-NS inspirals with a high accuracy, 
one always calibrates the model waveform by comparing 
it with the NR result. This implies that a high- accuracy 
numerical waveform is necessary. To achieve a high ac- 
curacy, one needs to perform more accurate simulations 
than the present ones. One of the keys for improving 
the accuracy is to reduce the orbital eccentricity in the 
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FIG. 9: The same as Fig.EJbut for MSI. 
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FIG. 8: The same as Fig.[7]but for H4. 



future. 
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FIG. 10: Snapshots of the orbital-plane density profile of NS-NS binaries in close orbits for APR4. The color denotes the 
density in units of logp(g/cm^). 




FIG. 12: Difference between (/!>nr and ^pn/eob for APR4-1215 (left panel) and for H4-1215 (right panel). Here we choose 
n = 1.8. 



